//--------------------------------------------------------------------------------------
//	 AceFile.hh
//---------------------------------------------------------------------------------------
// \verbatim 
// <pre>

//
// A note on coordinate systems.  There are three coordinate systems you must be 
// aware of to use this code.   These three are illustrated in the assembly 
// snippet below:
//
//            -             1         2         3            4
//            432101234567890123456789012345678901233334567890    Consed coords (depadded) 
// consensus       acttgtagaaagacacattccaaagatactaac***accttta...
// 
// r2         acttgtatttttgaaagacacattccaaagatactaacgggaccttta...
// r1              acttgtagaaagacacattccaaagatactaac***accttta... 
// r3                  gtagaaagacacattccaaagatactaac***accttta...
// r2a                                                  ccttta...
//            432101234567890123456789012345678901234567890123... Padded 1-based
//
//            543210123456789012345678901234567890123456789012... Padded 0-based
//            -              1         2         3         4 
//                           
// If you click on the first base of r2a in consed, the coordinate 
// reported will be the consed coordinate of 35.  If you look at the 
// AF line, however, the coordinate will be 38.  That is, the visual
// coordinates reported are Depadded coordinates.  The coordinates 
// stored in the AceFile (KJD: at least the AF lines, need to overtly 
// verify other coordinates) are Padded coordinates.  If you don't 
// realize this difference, unending confusion can follow.  The only time
// you need to be aware of the depadded coordinates are when you are 
// looking at the assembly in consed (KJD: actually, crossmatch results
// are given in depadded coordinates, so that is one other place where
// depadded coordinates arise). 
// 
// Everywhere in my code I try to use the padded coordinates.  However, 
// in my code there is the additional complication that I store everything
// in 0-based arrays.  Thus, the first consensus base has coordinate 0, 
// rather than 1.  The functions CIdx2RIdx and RIdx2CIdx convert from 
// AF coordinates found in the Ace file to my internal coordinates.  
//
// For example, for r2a, if I want the consensus position, in my code, 
// for base 0, Call RIdx2CIdx(r2a,0) that would yeild:  0+38-1 = 37.  
// 
// This is such an unending source of confusion that I am tempted to 
// switch my code to 1-base just to reduce the number of different 
// coordinate systems by one. 
// 
// -------------------
// AF and mAlignStart offsets for the above assembly:
//
// r2     -4
// r1      1
// r3      5
// r2a     38 
//
// </pre>
// \endverbatim
//
//
//---------------------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/include/AceFile.hh#3 $
//---------------------------------------------------------------------------------------
// $Log:$
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//

#ifndef _ACE_FILE
#define _ACE_FILE

#include "LocalSTL.hh"
#include "Sequence.hh"
#include "HashUtils.hh"
#include "FASTAFile.hh"
#include "PrintUtils.hh"
#include <iostream>
#include <iomanip>

class BaseSegment{
 public:
  int mStart;
  int mEnd;
  string mName;
};

//------------------------------------------------------------------------------
// * TemplateInfo
//------------------------------------------------------------------------------
// Class to hold information about template pairs.   
// 
class TemplateInfo{
 public:

  TemplateInfo(){
    mbOrientationOK = false;
  }

  string mTID;
  int    mLength;
  int    mInsertSize;

  bool   mbOrientationOK;
  bool   mbPair;
  bool   mbOverlapping;   // True if pair actually overlaps. 
  
  string mRName1;
  string mTDir1;
  char   mRDir1;
  int    mRSize1;
  int    mRStart1;

  string mRName2;
  string mTDir2;
  char   mRDir2;
  int    mRSize2;
  int    mRStart2;

  void Print(ostream &fout){
    fout<<mTID<<" ("<<mLength<<") "<<mbOrientationOK<<" "<<mbPair;
    fout<<"\t"<<mRName1<<" ("<<mTDir1<<","<<mRDir1<<")\t"<<mRStart1<<"\t"<<mRSize1;

    if (mbPair){
      fout<<"\t"<<mRName2<<" ("<<mTDir2<<","<<mRDir2<<")\t"<<mRStart2<<"\t"<<mRSize2;
    }
    fout<<endl;
  }
};

//------------------------------------------------------------------------------
// * Contig
//------------------------------------------------------------------------------
// Class representing a single contig within an Ace file.   Includes 
// the list of reads, the reads quality, trimming information, and a consensus.  
// 
class Contig : public Sequence{
 public:

  Contig(){mPrintLineSize=100;};

  int size(){
    return(Sequence::size());
  }
  int GetNumReads(){
    return(mReads.size());
  }

  void CalculateTrim(int cutoff);
  void CopyContig(Contig &C);

  void GetTemplateInfo(vector<TemplateInfo> &TInfo);
  void RebuildTemplateInfo();

  bool TemplateCoverIsZero();

  // Compute various coverage levels. 
  void GetTemplateCoverage(vector<int> &TCover);
  void GetBadTemplateCoverage(vector<int> &TCover);
  void GetHQDescrepency(vector<int> &HQDesc);
  void GetTemplateSpan(vector<int> &TCover);
  //  void Get10KTemplateCoverage(vector<int> &Cover);
  //void Get2KTemplateCoverage(vector<int> &Cover);
  void GetCoverage(vector<int> &Coverage);
  void GetWGSBACCoverage(vector<int> &WGSCover,vector<int> &BACCover);


  void ChooseBestTemplatePair(vector<int> &ReadIdxs,int &spanOut);
  void ComputeSyntheticBS();
  void PrintTemplateInfo(ostream &fout);
  void RecordTemplateInfo(int ReadIdx);
  void BuildTemplateStartEnd(vector<TemplateInfo> &TInfo,
			     Int2IntMap &startmap,Int2IntMap &endmap);

  void BuildBadTemplateStartEnd(vector<TemplateInfo> &TInfo,
				Int2IntMap &startmap,Int2IntMap &endmap);

  void BuildTemplateSpanStartEnd(vector<TemplateInfo> &TInfo,
				 Int2IntMap &startmap,Int2IntMap &endmap);

  void PrintAlignment(ostream &fout,int start,int end);
  void PrintAlignmentWithQuality(ostream &fout,int start,int end);

  int PrintBACTrimmedContig(ostream &fout,int minLength);

  int CalcCoverage(int idx);
  int CalcP20Coverage(int idx);
  int CalcCoverageTrimmed(int idx);

  long CountP20Bases();

  float CalcAverageTrimmedCoverage();
  float CalcAverageCoverage();
  float CalcAverageP20Coverage();

  int GetMatchingBasesTrimmed(int readIdx); 
  int GetMatchingBases(int readIdx); 

  int GetNumReadsInWindow(int winStart,int winEnd);
  string GetWindowReadName(int winStart,int winEnd,int rIdx);

  char GetAlignedCharacter(int readNum,int ConsensusPosition);
  string GetAlignedCharacters(int winStart,int winEnd,int ConsensusPosition);

  void GetTrimmedDepadded(Sequence &Seq);
  int    mPrintLineSize;

  void WriteAce(ostream &fout);

  void ComputeConsensus(int sIdx,int eIdx);
  void RemoveTails(int leftTail,int rightTail,vector<string> &RemovedReads);
  bool HasMateInContig(string R1);
  bool HasGoodMateInContig(string R1);
  void RemoveRead(string readName);


  void CopyConsensus(int targetStart,Contig &C,int sourceStart,int sourceEnd);
  void ComputeConsensus(int cIdx,char &consensus,int &quality);
  void ComputeConsensus2(int cIdx,char &consensus,int &quality);
  void CreateConsensusForCopy(); // Deprecated. 
  void FixupAlignStarts(int newStart,int &oldLeftMost,int &oldRightMost);
  
  int  CountReadsMatching(string key);
  bool AllReadsMatch(string key);
  bool NoReadsMatch(string key);

  void TrimPadTail();

  void ReverseComplementContig();

  void FindConsensusChangePoints(vector<string> &readList,int &leftmost,int &rightmost);

  void FindReadSets(int leftCrash, int rightCrash, 
		    vector<string> &LeftVec,
		    vector<string> &MiddleVec, 
		    vector<string> &RightVec);  
  
  void SplitContigRight(vector<string> &KillVec,
			vector<string> &RightVec,
			int leftmostAlteredConsensus,
			int rightmostAlteredConsensus,Contig &NewC);
  
  void SplitContigRight(vector<string> &RightVec,								 
			Contig &NewC);
  
  void MergeContigRight(Contig &C,int StartLocation);
  
  void SplitContig(int leftCrash, int rightCrash,Contig &NewC1,Contig &NewC2);
  
  void CopyRead(Contig &NewC,string ReadName);

  void AddReads(SequenceVector &Reads,vector<int> &offsets);
  void AddReadsNoBS(SequenceVector &Reads,vector<int> &offsets);
  void AddReadsNoBS(SequenceVector &Reads,vector<int> &offsets,vector<char> &senses);

  void  AddReads(SequenceVector &Reads,SequenceQualityVector &ReadsQual,vector<int> &offsets);


  // KJD Is AF in padded or depadded coordinates?  Check....
  
  //----------------------------------------
  // * CIdx2RIdx
  //----------------------------------------
  // \verbatim 
  // <pre>	
  //
  //            -             1         2         3            4
  //            432101234567890123456789012345678901233334567890    Consed coords (depadded) 
  // consensus       acttgtagaaagacacattccaaagatactaac***accttta...
  // 
  // r2         acttgtatttttgaaagacacattccaaagatactaacgggaccttta...
  // r1              acttgtagaaagacacattccaaagatactaac***accttta... 
  // r3                  gtagaaagacacattccaaagatactaac***accttta...
  // r2a                                                  ccttta...
  //            432101234567890123456789012345678901234567890123... Padded 1-based
  //
  //            543210123456789012345678901234567890123456789012... Padded 0-based
  //            -              1         2         3         4 
  // 
  // AF and mAlignStart offsets for the above assembly:
  //
  // r2     -4
  // r1      1
  // r3      5
  // r2a     38 
  //
  // So, to convert MY ReadIdx=0 to the read Idx for r2, 
  // 
  // AF[r2] = -4 -> CIdx = -5.  
  // AF[r1] = 1 -> CIdx = 0
  // AF[r2a] = 38 -> CIdx = 37. 
  //
  // AF[r2][20] ==> CIdx = 20 - 4 -1 = 15. 
  // 
  // </pre>
  // \endverbatim 
  //
  inline int RIdx2CIdx(int rIdx,int readNum){
    return(rIdx + mAlignStart[readNum]-1);
  }

  //----------------------------------------
  // * RIdx2CIdx
  //----------------------------------------
  // I use zero based coordinates.  Phrap uses one based 
  // coordinates.  So, as part of interpreting AlignStart values
  // this method also translates between coordinates. 
  // 
  //                        A G G G T A ...
  // MyContigIdx   -3 -2 -1 0 1 2 3 4 5 ...
  // 
  // PhrapCIdx     -2 -1  0 1 2 3 4 5 6 ...
  // 
  //                      N A G G G T A ...   mAlignStart = 0
  //                N  N  N A G G G T A ...   mAlignStart = -2
  //                          G G G T A ...   mAlignStart = 2
  // 
  // contigIdx is here taken to be MyContigIdx, which is zero-based. 
  // readIdx is here taken to be a zero-based readIdx. 
  // 
  // rIdx = 4.  mAlignStart = 0  cIdx = 3
  // rIdx = 6   mAlignStart = -2 cIdx = 3
  // rIdx = 2   mAlignStart = 2  cIdx = 3


  inline int CIdx2RIdx(int contigIdx,int readNum){	 
    return(contigIdx - mAlignStart[readNum]+1);
  }


  // private: 

  //string mName;
  int mPaddedBases;
  int mNumEntries;
  int mNumBaseSeqs;
  
  String2IntMultiMap   mTemplate2ReadMap;   // Maps template names to associated reads. 
  String2StringMap     mRead2TemplateMap;

  // Maps read idx to template dirs.  I made this a map rather than a vector 
  // since the way I parse the template info make it a slight pain to ensure
  // that there is one entry per read. 
  Int2IntMap   mReadTemplateDir; 

  // Returns the index of a read based on the name. 
  String2IntMap    mReadIdxMap;

  //Sequence              mConsensus;
  SequenceQuality       mConsensusQuality;

  SequenceVector        mReads;
  SequenceQualityVector mReadsQuality;  

  vector<TemplateInfo> mTInfo;

  vector<char>          mTemplateDir;

  // Auxillary 
  vector<int>           mAlignStart;
  vector<char>          mReadSense;
  vector<int>           mReadTrimLeft;
  vector<int>           mReadTrimRight;

  vector<BaseSegment>   mBaseSegments;  
  vector<string>        mQALines; // KJD Unparsed for now. 

  int                   mTrimLeft;
  int                   mTrimRight;

  int    mLeftmostBAC;
  int    mRightmostBAC;

  // The trim which is the most interior of the first BAC or the 
  // quality trim. 
  int    mBACTrimLeft;
  int    mBACTrimRight;

  // A bit of a hack.  When a new contig is created by split, this stores 
  // the leftmost and rightmost contig position of any read in the 
  // original contig.   
  int    mOldLeftMost;
  int    mOldRightMost;

};

//------------------------------------------------------------------------------
// * AceFile
//------------------------------------------------------------------------------
// Class to almost fully parse a new-ace format Ace file. 
// 
// Some information about the format of the ACE file gleened by examining the 
// phrap code is given here for general reference:
//
// AS (num_contigs-(numsingletons+num_duplicates)  num_reads
// CO  ContigName padded_bases contig->num_entries num_base_seqs
// {contig sequence}
// AF  ReadName  reverse?->C/U  align_entry->start 
// AF  ReadName  
// ...
// BS  start end  ReadName
// BS  start end  ReadName
// RD  ReadName padded_bases whole_read_infos read_tags
// {read sequence}
// QA  l_pad_qual_trim  r_pad_qual_trim  l_pad_align_trim+1 r_pad_align_trim+1
// DS  description line.  
// RT  id(align_entry->seq_entry)  tag->type  -1-tag->start  -1-tag->end params->date
//
// Some additional Ace file facts:
//
// Consensus bases for a contig are 1 based in the file. 
// Consensus coordinates in Consed are shown as if without pads. 
// 
// The RD lines in an Ace file are already in the correct sense for the alignment
// (as opposed to the original sequence in the fasta file which are in their 
// original sense).  
//
// The RD lines include pads needed for alignment. 
//
// KJD: What are WR lines:
//
// WR{
// 2kContig140_Contig33R1U0000106 unaligned phrap 030218:142945
// }
// 
// Does my code choke when it sees this?  I only recently encountered such lines. 
// 
//
//
// Facts about my code:
//
// I use zero based sequences every where.  The first consensus base is zero. 
// The first base of a read is zero, etc.   I convert to Ace coordinates as needed. 
// Mainly this conversion occurs whenver I use the AF values, which I preserve 
// unaltered in mAlignStart, rather than translated.  It'd probably make more sense
// to translate the AF lines as I read/write them, then ALL my own code would be 
// zero based and the translation would only occur on read/write, but for now 
// the AF lines are kept in original coordinates and translated in CIdx2RIdx/RIdx2CIdx
//
// The AceFile class is simply a list of Contigs with some overall stats and helper
// functions.  Each contig is a consensus, a consensus quality, a list of reads, 
// optionally a list of read qualities, and some convenience data structures to make
// it easier to navigate.  Be careful that the convenience data structures are 
// initialized before you try to use them!  Some are initialized on Read(), others
// require explicit initialization. 
// 
class AceFile : public vector<Contig>{

 private:
  string mName;
  SequenceQualityVector  mContigQualityVec;

  // This is a multi-map because, unfortunately, reads can occurr in more 
  // than one contig.  In general they shouldn't, of course, but in practice 
  // they sometimees do and we need to handle that case rationally. 
  String2IntMultiMap    mReadContigMap;

  bool   mbIsConsistent;
  int    mNumReads;

  // Stats
  int    mNumContigs;                        
  int    mTotalReadCount; 
  int    mWGSReadCount;                         // Whole Genome Shotgun Read count. 
  int    mContigsMeetingCriteria;
  int    mTrimmedContigsMeetingCriteria;
  int    mSizeMeetingCriteria;
  int    mTrimmedSizeMeetingCriteria;
  int    mN50;
  int    mN50Trimmed;
  float  mAverageCoverage;
  float  mTrimmedAverageCoverage;
  int    mTotalContigLength;
  int    mTrimmedContigLength;
  int    mMaxContigSize;
  int    mMaxTrimmedContigSize;

  // Parameters
  //  int    mPhredTrim;
  //int    mContigTrim;
  int    mMinContigLength;

  // A map that tells which contig contains reads with given TID
  //  String2IntMultiMap mTIDContigMap;  

  //====================

 public:

  int    mPhredTrim;
  int    mContigTrim;
  
  // This is sugar for Java for right now.  Eventually, I will want to 
  // move the vector<Contig> part of AceFile inside a simple class and 
  // overload operator [] to get backward compatibility. 
  int  size(){
    return(vector<Contig>::size());
  }

  Contig& GetContig(int c){
    return((*this)[c]);
  }


  bool gzExists(string FileName); // convenience method to detect if zipped file exists. 

  void Read(istream &fin);
  void Read(string AceName);
  void ReadAll(string &RootName);
  void ReadAll(string &AceName,string &ReadQualName);

  int  CName2CIdx(string ContigName){
    for(int c = 0;c< size();c++){
      if (Contains((*this)[c].mName,ContigName)) return(c);
    }
    return(-1);
  }
  
  void Write(string FileName){
    ofstream fout;
    fout.open(FileName.c_str());
    int numReads = 0;
    for(int i = 0;i < size();i++){
      numReads+=(*this)[i].mReads.size();		
    }

    fout<<"AS "<<size()<<" "<<numReads<<endl;
    for(int c = 0;c< size();c++){
      (*this)[c].WriteAce(fout);
    }
    fout.close();
  }

  void Write(string FileName,StringSet &SkipList){
    ofstream fout;
    fout.open(FileName.c_str());
    int numReads = 0,numContigs=0;
    for(int i = 0;i < size();i++){
      if (!SkipList.Contains((*this)[i].mName)){
	numContigs++;
	numReads+=(*this)[i].mReads.size();		
      }
    }

    fout<<"AS "<<numContigs<<" "<<numReads<<endl;
    for(int c = 0;c< size();c++){
      if (!SkipList.Contains((*this)[c].mName)){
	(*this)[c].WriteAce(fout);
      }
    }
    fout.close();
  }

  void WriteAceFromContigs(string OutFile,vector<Contig> &OutputContigs){
    string ContigsFile = OutFile+".contigs";
    FASTAFile FA;
    
    ofstream fout,ctgout;
    fout.open(OutFile.c_str());
    if (fout.fail()){
      cout<<"Unable to open "<<OutFile<<endl;
    }
    ctgout.open(ContigsFile.c_str());
    
    int totalContigs = OutputContigs.size();
    cout<<totalContigs<<" being written."<<endl;
    // figure out total reads. 
    int totalReads = 0;
    for(int i = 0;i< OutputContigs.size();i++){
      totalReads +=OutputContigs[i].mReads.size();
    }    
    fout<<"AS "<<totalContigs<<" "<<totalReads<<endl;
    int globalContigNum = 0;
    for(int i = 0;i< OutputContigs.size();i++){
      string CName = "Contig"+itoa(globalContigNum++);
      OutputContigs[i].mName = CName;
      OutputContigs[i].WriteAce(fout);
      
      FA.Write(ctgout,OutputContigs[i]);
    }    
    fout.close();
  }


  void WriteGZipped(string FileName){
    ozipstream fout;
    fout.open(FileName.c_str());
    int numReads = 0;
    for(int i = 0;i < size();i++){
      numReads+=(*this)[i].mReads.size();		
    }

    fout<<"AS "<<size()<<" "<<numReads<<endl;
    for(int c = 0;c< size();c++){
      (*this)[c].WriteAce(fout);
    }
    fout.close();
  }

  bool FindNewestAce(string &AceRoot);
  void ComputeAllContigStats();
  void CleanInitFromContigs();
  void CalcAverageCoverage();
  void CalcAverageCoverageTrimmed();
  void AddFakeQualityValues();
  void CopyBQToQuality();
  void PrefixContigs(string prefix);

  // These are some functions to aid in accumulating statistics across 
  // several 
  // Sorry this is here.  ContigQualVec is here, so that's where this needs to be, for now. 
  //  void PrintAlignmentCoverageTrimmedWithQual(int contigNum,ostream &fout,int start,int end);


  void AddReadsQuality(SequenceQualityVector &Quals);
  void Read(string FName,SequenceQualityVector &ReadQuals);
  void GetReadsQuality(istream &fin);
  void BuildReadHash();

  void CalculateReadTrims(int cutoff);
  void CalculateContigTrims(int cutoff);
  void CalculateBACContigTrims(int cutoff);


  int  N50();
  void GetSortedSizes(vector<int> &sizeList);

  int  N50Trimmed();
  void GetSortedSizesTrimmed(vector<int> &sizeList);

  bool GetNextQuality(istream &fin,SequenceQuality &Qual);

  bool IsBACName(string Name);


  void RemoveContig(int c);

  void FindTemplatesBetween(Contig &C1,Contig &C2,vector<TemplateInfo> &TInfo);
  void FindTemplatesBetween(int c1,int c2,vector<TemplateInfo> &TInfo);
  void GetContigsSharingTemplates(int cIdx,vector<int> &ConnectedContigs);

  void MultiSplitContig(int MinMiddle, int firstCIdx, vector<int> &leftSplits,vector<int> &rightSplits);
						
  void AdjustSplitPoints(int currentPair, int firstCIdx,int newCIdx,int newCIdx2,
			 vector<int> &leftSplits,vector<int> &rightSplits,
			 vector<int> &SplitContig);
  
 // ------------------ JAVA INTERFACE SUGAR -----------------------
  // 
  // Routines to make it easier to show progress reading in a file from 
  // Java, and also a way to isolate izipstream from the Java SWIG interface.
  // 
  izipstream mZin;
  int mContigNum;
  string mLine;
  int Open(string FName);
  void Close();
  int ReadContigBatch();
  // ----------------------------------------------------------------
 
 protected:

  // Parse functions for the individual sections of the ACE File. 
  string ParseCO(istream &fin,string &line,int &ContigNum);
  string ParseBQ(istream &fin,string &line,int &ContigNum);
  string ParseAF(istream &fin,string &line,int &ContigNum);
  string ParseBS(istream &fin,string &line,int &ContigNum);
  string ParseRD(istream &fin,string &line,int &ContigNum);

  /*
    String2IntMap    mTemplateMap;
    vector<TemplateInfo>  mTemplateInfo;
    void RecordTemplateInfo(string &name,string &desc,int contigIdx,int alignStart,
    int readSize,char sense);
  */
};

#endif 


